NBI-HE-00-24 
NORDITA-2000/50-HE 
hep-th/00051^ 
May 2000 



Monte Carlo Studies of the I IB Matrix Model at Large N 



J. AmBJ0RN^)[], K.N. ANAGNOSTOPOULOS^^fl, 
W. BlETENHOLZ^)0, T. HOTTA^)0 AND J. NlSHIMURA^^Q 

Niels Bohr Institute, Copenhagen University, 
Blegdamsvej 17, DK-2100 Copenhagen 0, Denmark 
^) Department of Physics, University of Crete, 
P.O. Box 2208, GR-71003 Heraklion, Creece 
3) NORDITA 
Blegdamsvej 17, DK-2100 Copenhagen 0, Denmark 
^) Institute of Physics, University of Tokyo, 
Komaha, Meguro-ku, Tokyo 153-8902, Japan 



Abstract 

The low-energy effective tfieory of the IIB matrix model developed by H. Aoki et al. 
is written down explicitly in terms of bosonic variables only. The effective theory is 
then studied by Monte Carlo simulations in order to investigate the possibility of a 
spontaneous breakdown of Lorentz invariance. The imaginary part of the effective 
action, which causes the so-called sign problem in the simulation, is dropped by hand. 
The extent of the eigenvalue distribution of the bosonic matrices shows a power-law 
large behavior, consistent with a simple branched-polymer prediction. We observe, 
however, that the eigenvalue distribution becomes more and more isotropic in the 
ten-dimensional space-time as we increase A^. This suggests that if the spontaneous 
breakdown of Lorentz invariance really occurs in the IIB matrix model, a crucial role 
must be played by the imaginary part of the effective action. 
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1 Introduction 



For more than two decades superstring theories have been studied as the most promising 
candidates for a unified theory of all the interactions including gravity. The theories have 
the potential to predict the space-time dimensionality, the gauge group, the matter content, 
and so on, from first principles. The existence of infinitely many perturbative vacua implies, 
however, that an understanding of nonperturbative effects is crucial to extract information 
about the real vacuum of the theory. Recent proposals for nonperturbative formulations of 
superstring theories [l], |, H, |^ may therefore be of analogous importance for our under- 
standing of non-perturbative aspects of string theory as lattice gauge theory has been 
in understanding nonperturbative dynamics of gauge theories. The IIB matrix model 0, 
which is conjectured to be a nonperturbative formulation of type IIB superstring theory 
(for a review, see Ref. |]^), takes the form of a large reduced model and in the same 
way that Monte Carlo studies of lattice gauge theory clarified many important nonpertur- 
bative aspects of the strong interaction, Monte Carlo studies of the IIB matrix model might 
illuminate the nonperturbative dynamics of superstring theories^. A number of numerical 



studies have already been done [|T^, [11], |T2|' HI; HI' 0' 0) O HI' Hll pursue that direction. 
Earlier numerical studies of world-sheet perturbative aspects of superstring-like theories can 
be found in Refs. Q ^ 0]. 

In this paper, we study the IIB matrix model at large by Monte Carlo simulations. In 
particular, we investigate the possibility of a spontaneous symmetry breakdown of Lorentz in- 



variance using the low-energy effective theory of the IIB matrix model proposed by Ref. [23 



The authors of also study this issue, using an approach where the bosonic and fermionic 
matrices are both decomposed into diagonal elements and off-diagonal elements, and the off- 
diagonal elements are integrated out first perturbatively. Such a perturbative expansion is 
valid when the bosonic diagonal elements, which can be regarded as coordinates of A^ points 
in ten-dimensional fiat space-time, are well separated from each other. In other words, what 
one obtains after integrating over the off-diagonal elements perturbatively can be considered 
as a low-energy effective theory of the IIB matrix model. Note that this perturbative expan- 
sion has nothing to do with the string perturbative expansion with respect to the world-sheet 
topology. Therefore, even at the one-loop level, the low-energy effective theory is expected 
to include nonperturbative effects of superstring theory, provided, of course, that the IIB 



^ Monte Carlo studies of Matrix Theory |^] would be technically more involved due to a lattice discretiza- 
tion of the time direction iSl. 
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matrix model conjecture is true. 

In fact, in order to obtain the low-energy effective action only for the bosonic diagonal 
elements, one still has to perform the integration over the fermionic diagonal elements, which 
is nontrivial, since their action turns out to be quartic. Although the explicit form of the 



final low-energy effective theory has not been derived, the theory was shown to be 
described by some complicated branched-polymer like system, typically involving a "double- 
tree" structure, in a flat ten-dimensional space time. Thus, even at the one-loop levelQ, the 
low-energy effective theory contains highly nontrivial dynamics. It was further argued that 
the double-tree structure of the one-loop low-energy effective theory might be responsible 
for a collapse of the distribution of the bosonic diagonal elements. The first Monte Carlo 
results of a branched-polymer system with a double-tree structure was reported in Ref. [0. 

Here, we write down explicitly the low-energy effective theory of the IIB matrix model in 
terms of bosonic variables only. Instead of integrating over both bosonic and fermionic off- 
diagonal elements first, we leave the bosonic off-diagonal elements unintegrated. The action 
for the fermionic diagonal elements is then still quadratic and can be integrated explicitly, 
yielding a Pfaffian. Integration over the bosonic off-diagonal elements as well as the bosonic 
diagonal elements can be done by Monte Carlo simulation. In other words, the bosonic 
off-diagonal elements play the role of auxiliary variables, which enable us to simulate the 
complicated branched-polymer like system describing the dynamics of the bosonic diagonal 
elements. 

The Pfaffian induced by the integration over the fermionic diagonal elements is generically 
complex. In general, when the action of a theory has a non-zero imaginary part, the number 
of configurations needed to extract any information increases exponentially with the system 



size, except in a few situations, where alternative sampling methods can be invented p^ , p5[| . 
This notorious technical problem in Monte Carlo simulations is known as the "sign problem" . 
In fact, the problem exists already in Monte Carlo simulations of the original IIB matrix 
model and it is inherited by the low-energy effective theory. In the present work, we simply 
use the absolute value of the Pfaffian in order to avoid the sign problem and examine only 
the effect of the modulus of the Pfaffian. Our results suggest that there is no spontaneous 
symmetry breakdown (SSB) of Lorentz invariance. From this we conclude that if the SSB 
ever occurs in the IIB matrix model, the phase of the Pfaffian must play a crucial role. 



^The validity of the one-loop approximation for studying the low-energy dynamics of supersymmetric 
large TV reduced models has been demonstrated in Ref. ||l8| through the study of the four-dimensional 
version of the IIB matrix model. We will come back to this point is Section ^. 
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We also study the six-dimensional version of the JIB matrix model for comparison. The 
conclusion is the same, but we find an intriguing difference in the finite effects. 

The paper is organized as follows. In Section ^, we describe the definition of the IIB 
matrix model and review some important properties relevant for this work. In Section ^ 
we derive the low-energy effective theory of the IIB matrix model and explain the model 
we investigate by Monte Carlo simulations. In Section ^, we present our results for the 
distribution of the bosonic diagonal elements. In particular, we discuss the possibility of SSB 
of Lorentz invariance. Section ^ is devoted to a summary and conclusions. In Appendix A, we 
explain the details of the algorithm we use for the Monte Carlo simulation. In Appendix B, 
we present some systematic studies for optimizing the parameters involved in the algorithm. 



2 The IIB matrix model 



The IIB matrix model is formally a zero-volume limit of ten-dimensional pure M = 1 
supersymmetric Yang-Mills theory. The action, therefore, is given by 
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(yU = 1, ■ ■ ■ , 10) and ^/'q, (a = 1, ■ ■ ■ , 16) are N x N traceless Hermitian matrices, which 
can be expanded in terms of the generators t° of SU(A^) as 



(2.4) 



a=l 
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where is a real variable and ip'^ is a real Grassmann variable. We assume that the 
generators are normalized as tr(t"t'') = 6ab- The measure dip in ( p.lD is defined by 
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and similarly for dA. The model (|2.1|) appears after a Wick rotation, so that the metric has 



Euclidean signature. The 16 x 16 matrices are defined by 

f^=Cr^, (2.6) 

where are ten-dimensional gamma matrices after Weyl projection, and the unitary matrix 
C is a charge conjugation matrix satisfying 

cr^ct = (r^)^ ; c^ = c. (2.7) 

Due to (|2.7| ), the matrices are symmetric. 

An explicit representation of the gamma matrices is given by 

Fi = i (72 ® (72 ® CT2 ® 0-2 ; F2 = i 0-2 (g) (T2 ® 1 ® 0^1 ; F3 = i (T2 ® 0-2 ® 1 ® CTs ; 

F4 = z 0-2 ® 0-1 (g) 0-2 (8) 1 ; F5 = z 0-2 ® 0-3 (g) 0-2 ® 1 ; Fg = i 0-2 (8 1 (g ai 0-2 ; 

F7 = 2o-2®l®cr3(g)a2; F8 = io-i(g)l(g)l(g)l; Tg = i ® 1 ® 1 ; 

Fio = l®lOl®l, (2.8) 

for which the charge conjugation matrix C becomes a unit matrix and therefore F^ = F^. 

The model has a manifest ten-dimensional Lorentz invariance, by which we actually mean 
an SO (10) invariance. transforms as a vector and ipa transforms as a Majorana-Weyl 
spinor. The model is manifestly super symmetric, and also has a SU(A^) symmetry 

A,^VA,V^ ■ ^^^V^^V^ , (2.9) 

where V G SU(A^). All these symmetries are inherited from the super Yang-Mills theory 



before taking the zero- volume limit. In particular, the SU(A^) symmetry (|2.9|) is a remnant 
of the local gauge symmetry. 

The fermion integral Zi[A\ in (|2.1| ) can be obtained explicitly as 



Zi[A] = FiM , (2.10) 



where 
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fabc = -^tr(r[^^^^]) . (2.11) 



4 



The real totally-antisymmetric tensor fabc gives the structure constants of SU(A^) and the 
matrix M.aa,bi3 is a 16 (A^^ — 1) x 16 (A^^ — 1) anti-symmetric matrix, regarding each of {aa) 
and {h(3) as a single index. 

The convergence of the integration over the bosonic matrices in ( |2.1|) is nontrivial since 
the integration domain for the Hermitian matrices is non-compact. Even for finite there 
is a potential danger of divergence when the eigenvalues of become large. This issue 
has been addressed in Ref. |]23[ using one-loop perturbative arguments which pointed to the 
finiteness of the IIB matrix model given by (|2.1j ) for arbitrary A^. This conclusion is in 
agreement with an exact result available for N = 2 ||3l| and a numerical result obtained 



for A^ = 3 [|T2|. Thus it is conceivable that the above conclusion, obtained by one-loop 
arguments, holds in general. 

Since the model is well-defined without any cutoff, the parameter g, which is the only 
parameter of the model, can be absorbed by rescaling the variables, 

= g"^X,, (2.12) 
= . (2.13) 

Therefore, g is & scale parameter rather than a coupling constant, i.e. the g dependence of 
physical quantities is completely determined on dimensional grounds^. In what follows, we 
take g = 1 without loss of generality. 

For comparison, we also study the six-dimensional version of the IIB matrix model. In 
this regard we recall that pure A/" = 1 supersymmetric Yang-Mills theory can be also defined 
in 3D, 4D and 6D, as well as in lOD. Hence, by taking a zero-volume limit of these theories, 
we arrive at supersymmetric large A^ reduced models which are D = 3,4,6 versions of the IIB 
matrix model. Using the one-loop argument mentioned above, one concludes that the model 
is ill-defined for D = 3, but well-defined for D = 4,6, 10, irrespectively of A^. For D = 4, 
Monte Carlo simulations up to A^ = 48 further confirms this statement [|l8l. The effective 
action induced by fermions (logarithm of the fermion integral ^f[^]) is real for D = 4. It is 
complex in general for D = 6 and D = 10, however, which causes the sign problem in Monte 
Carlo simulations. 



•^The scale parameter g should be tuned appropriately as one sends N to infinity so that all the correlation 
functions of Wilson loops have a finite large N limit. Whether such a limit really exists or not is one of 
the important dynamical issues, which was addressed in Ref. [p^ for the four-dimensional version of the IIB 
matrix model. 
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3 Low-energy effective theory 



In this section, we derive the low-energy effective theory of the IIB matrix model (|2.1| ) along 



the lines discussed in Ref. |£3|. We first decompose the N x N Hermitian matrices and 

ipa as 

i)a = ia + ^a , (3.1) 

where and represent the diagonal parts, while and (fa represent the off-diagonal 
parts. We also introduce N ten-dimensional vectors Xi {i = 1, ■ ■ ■ , N) through = 
and N ten-dimensional Majorana-Weyl spinors {i = 1, - ■ ■ ,N) through = {ia)ii- For 
the off-diagonal elements, we use the notations a^j-,- = (a.^)ij and (paij = {0a)ij, where i ^ j. 
Using the decomposition the actions ( |2.2|) and ( p.3|) can be written as 



Sh = -ti (^-^a^[x^,[Xf,,d^]]-^[Xf,,afJ^^ +[x^,a^][a , (3.2) 

5'f = (f^)a/3tr(^(^„[x^,<^^] - (^c.[^/3,aM] - + "^"[^M'*^/?]) • (3-3) 

The one-loop approximation amounts to keeping only the quadratic terms in a and (f in 
the above expressions, neglecting the higher-order terms, i.e., the O(a^) term and the O(a^) 
term in and the 0{aip'^) term in 5'f. The quadratic term in a in ( p.2| ) has zero modes 



due to the fact that the original model (p.l|) has the SU(A^) invariance ( p.9|) . We thus have 



to "fix the gauge" properly. Following Ref. [^], we choose the "gauge- fixing" term (which 



is the reduced model counterpart of the Feynman gauge in ordinary gauge theory) and the 
corresponding Faddeev-Popov ghost term as 

-^g.f. = -^tr([x^,a^]2) , 

5'gh = -tT{[x^,b][Xf, + a^,c]) . (3.4) 



in the ghost action can be neglected within the one-loop approximation. Integration 
over the ghost fields b, c then yields {A(x)}^, where A{x) is defined as 

A(x) = n(^^ - ^.)' • (3-5) 
Therefore, the partition function at the one-loop approximation can be written as 

7(l-l00p) 



^IIB 



J dxdae-^G{A(x)}2zf~'°°P)[x,a] , (3.6) 



where 

Sg = ^{xi - Xj^la^ijl'^ , (3.7) 

i<j 

and the one-loop approximated fermion integral is defined as 

zf-^°°P)[x,a] = y"dedy.e-^f'\ (3.8) 

5f ^ = {tf,)al3tT(^0a[x^, 0(3] " 0a[ip, Cl^] ' (^^] j . (3.9) 



In what follows we calculate zf [x, a] explicitly. We first note that Sf^ can be 
written as 



Sf^ = -\ ^0cji{ff,)ap{Xif, - Xj^)0pij - i "^^iaM^^jp^jp , (3.10) 

where (paij in (|3.1CI| ) are defined by 



[Xi Xj) 



and M^q,j73 is a 16A^ x 16A^ matrix given as 

ML,ji3 = ^^!^2 i^ ^^^'^^)a(s{at^jiaaij - a^ija^ji) for i^j, (3.12) 

[Xi Xj ) 

ML,p = -Y.ML,p- (3-13) 

Integration over (p can now be done yielding {A(x)}^. We then integrate out ^Na using the 
delta functions in (|2.5D, yielding a factor of 1/A^^ followed by a replacement 

CiVa^-J^Oa (3.14) 

in ( |3.1(J| ). The integration over the (i = 1, ■ ■ ■ , (A^ "I)) yields Pf M . The matrix M is a 
16 (A^ — 1) X 16(A^ — 1) complex matrix defined as 

M,o.,jp = M[^.p - M'^^.p - Mi^^^ + M'^^^^p , (3.15) 

where indices i and j run from 1 to — 1. Note that there are identities 

m;„,^ = ^«,,/3 , ^/3,,« = , (3.16) 
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and similarly for Miajp. This means in particular that Mia j 13 is an anti-symmetric matrix, 
regarding each of (ia) and {jf3) as a single index. Thus the one-loop approximated fermion 
integral is obtained as 

zt"''"\x,a] = {A{x)r^^FfM . (3.17) 
We have checked numerically that indeed 

Zi[A] ~ Zl^'^°"'^\x, a] ; (A^)ij = Xi^ Sij + a^ij , (3.18) 

holds when the Xj's are well-separated and a^jj are generated with the distribution e"'^'^, 
where Sq is given by ( p.7|) . Note that the size of the matrix M in ( p.l7|) is of O(A^), whereas 



the size of the matrix Ai in ( |2.1CI| ) is of 0(A^^). The huge reduction is essentially because 
the integration over the fermionic off-diagonal elements has been done explicitly yielding 
{A(x)}® in (|3.17| ). This is the substantial gain from using the one-loop approximation. 



First we note that the model ( p.6|) as it stands has a singularity]^ at a;j^ = due to 
the singularity in ( p.l2| ). Therefore if one simulates the model (|3.6|) , the distribution of 



Xi collapses to the origin. We recall that such an ultraviolet singularity does not exist in 
the original IIB matrix model. On the other hand, the one-loop approximation is valid only 
when X'l S 8bTQ well separated from each other. Therefore, we have to introduce a UV cutoff to 
the distribution of Xi in order to make the model meaningful. As a UV cutoff, we introduce 
in the action the term given as 



Scut = Y.f{^{x.-x,y) , (3.19) 

i<j 

where the function f{x) is taken to be 

The dimensionless "spring constant" k should be taken to be large enough to prevent the 
Xj's from coming closer to each other than the cutoff £ (See Figs. |I] and |^.). Thus we arrive 
at the low-energy effective theory of the IIB matrix model ( f^.lj ) 



Zleet = J dxdae-^'^-^^"' {A(a;)}^°PfM , (3.21) 



which is written in terms of bosonic variables only. We have omitted the irrelevant constant 
factor in ( p.l7|) . 



"^This can be seen more clearly by rescaling a^ij as in (A^). Then, one finds that all the Xi^^-dependence 



of the partition function is contained in the Pfaffian, which has the singularity. 
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In fact, the UV cutoff parameter I can be scaled away from the theory ( |3.21| ) by rescahng 
the variables This means that the dependence of the 

results on the UV cutoff i is determined completely on dimensional grounds. In particular, 
dimensionless quantities are independent of i. When we are interested in dimensionful 
quantities, and in particular in their large N behavior, we have to know the iV-dependence of 
the UV cutoff i. In Ref. |^, it was argued that i should be taken to be A^-independent, based 
on a reasonable assumption that the ultraviolet behavior of the space-time structure of the 
IIB matrix model is controlled by the SU(2) matrix model. The issue is also addressed in the 



four-dimensional version of the IIB matrix model in Ref. [|18|. There, an A^-independent UV 
cutoff was shown to be generated dynamically by treating the full model nonperturbatively 
instead of making perturbative expansions around diagonal matrices. We therefore take 
£ = 1 for all in the following. 

The Pfaffian Pf M in ( |3.21D is complex in general. This poses the notorious "sign prob- 
lem", when one tries to study the model by Monte Carlo simulation. Note, however, that 
the problem is simply inherited from the original IIB matrix model as can be seen 

from (|3.18 ). In the present work, we take the absolute value |PfM| and throw away the 



phase by hand. This corresponds to taking the absolute value |Pf of the Pfaffian in the 
IIB matrix model and studying its low-energy effective theory. 

To summarize, the model we simulate is given by the partition function 



Z = y"da;dae-^G-^="* {A(x)}^° |PfM| . (3.22) 



Sq and S'cut are given by (|3.7|) and ( p.l9|) respectively, and M is defined through (|3.12|) , (|3.13| ) 



and ( |3.15| ). Note that the model ( p.22| ) still has the lOD Lorentz invariance. We would like to 



investigate whether the lOD Lorentz invariance breaks down, say, to 4D Lorentz invariance. 
The model corresponding to the six-dimensional version of the IIB matrix model can be 
obtained similarly. 

The validity of the one-loop approximation for studying the low-energy dynamics of 
supersymmetric large N reduced models has been addressed in Ref. |jl8| by studying the four- 
dimensional version of the IIB matrix model without any approximations. It was found that 
the large N behavior of the extent of the distribution of Xi agrees with the prediction from the 
one-loop, low-energy effective theory. We therefore expect that a low-energy phenomenon 
such as SSB of Lorentz invariance can be studied with the low-energy effective theory. 

On the other hand we note that from a technical point of view the one-loop approximation 
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offers certain advantages which are essential in the present work. We recall that even in the 
four-dimensional version of the IIB matrix model, the largest one can achieve for the 
full model (without the one-loop approximation) is = 48 using supercomputers. In order 
to detect even a small trend of SSB of Lorentz invariance from lOD to 4D, say, one would 
expect that should be larger than 4^ = 256. Due to the fact that the size of the matrix 
M in ( |3.22 ) is order of smaller than the size of the matrix M. in ( 2.10 ), the one- loop 
approximation enables us to reach A^ = 512 with reasonable effort. Details of the algorithm 
used in the Monte Carlo simulations are presented in the Appendix A. The algorithm is a 
variant of the Hybrid Monte Carlo algorithm [^] which is one of the standard algorithms 
used in the study of systems with dynamical fermions. The computational effort of the 
algorithm is estimated to be 0(A^^)~0(A^^/^), which should be compared with an estimate 
0(A^^) for the algorithm used for the full model in Ref. |1H . 



4 Results for the distribution of Xi 

First we look at the distribution p(r) of the distance r between Xj's, where the distance 



between two arbitrary points Xi ^ xj is measured by a/ {xi — XjY- In Fig. 0, we plot the 
results for D = 10 with A^ = 192, 256, 384 and 512. Fig. | shows the results for D = 6 with 
A^ = 256, 512 and 768. We note that in both cases the distribution at small r falls off rapidly 
below r ~ £, where £ = 1 is the UV cutoff introduced in (|3.20|) . The small penetration into 



the r < I region is due to k being finite. However, the results show that the values of k 
we have taken (k = 300 for D = 10 and k = 100 for D = 6) are large enough to make the 
penetration reasonably small. 

In order to see the spontaneous breakdown of Lorentz invariance, we consider the moment 
of inertia tensor of A^ points Xi {i = 1, ■ ■ ■ , N) in a fiat D-dimensional space-timeQ . It can 
be defined as 

2 

T^u = nT( ]\T _ 1 \ y^(^»A' ~ ^jli){^iu — Xjy) , (4-1) 

which is a. D X D real symmetric matrix. The D eigenvalues Ai > A2 > ■ ■ • > Ad > of the 
matrix T represent the principal moments of inertia. We measure A^ for each configuration 
and take an average (A^) over all the configurations generated by the Monte Carlo simulation. 
^Such a quantity has also been studied in Refs. p3| . 



10 



0.12 p 

0.1 - 

0.08 - 

0.06 - 

0.04 - 

0.02 - 

- 



T 1 1 

N = 192 I — \ — I 
N = 256 
N = 384 

N = 512 ^ □ i 



0.5 




Figure 1: The distribution of distances p(r) for D = 10 is plotted against r for N 
256, 384 and 512. 



192, 



If Lorentz invariance is not spontaneously broken, (A^) must be all equal in the large limit, 
representing an isotropic distribution of Xj. We therefore search for a trend which differs 
from such a large N behavior. 

We first note that if the system that describes the dynamics of Xi were a simple branched 
polymer in a flat /^-dimensional space-time, spontaneous breakdown of Lorentz invariance 
would certainly not occur, and moreover, the large N behavior of (A^) could be expected to 
be (A^) ~ A^^/^. This is due to the fact that the Hausdorff dimension of a branched polymer 
is (in = 4 irrespectively of the dimension D of the space-time in which it is embedded. For 
this we note that the extent of the distribution of Xi is given by 



R 



\ 



D 



(4.2) 



which is related to the number of points through ~ {R/iY^. The i is the UV cutoff 
taken to be £ = 1. Therefore, we find R ~ A^^/'^h ^ N'^/^^ which leads to the announced 
large A^ behavior (A^) ~ N'^^'^. 

In Fig. ^ we plot the normalized eigenvalues (A^) /N^^'^ against 1/A^ for D = 10 with A^ = 
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Figure 2: The distribution of distances p(r) for = 6 is plotted against r for = 256, 512 
and 768. 

192,256,384,512. We find that the smallest (normalized) eigenvalue (Aio)/A^^''^ is almost 
constant and the larger ones are monotonously decreasing towards the same constant. Thus 
the observed large N behavior suggests that there is no spontaneous breakdown of Lorentz 
invariance and the large N behavior of the extent of the x-distribution is consistent with 
a simple branched-polymer prediction. In particular, we see no trend for a gap developing 
between the fourth and the fifth largest eigenvalues, which could have been observed if the 
lOD Lorentz invariance were broken down to a four- dimensional one. In Fig. ^ we present 
the results for D = 6 with N = 192, 256, 512, 768. The qualitative behavior is the same as 
inD = 10. 

The results for D = 10 might suggest that the leading finite effect is given by 1/A^. 
For example, a linear extrapolation in 1/A^ using the data for = 384, 512 leads to almost 
the same results for all the ten (normalized) eigenvalues at = oo. It is therefore tempting 
to speculate that finite A^ effects in the IIB matrix model is given by a 1/A^ expansion. Such 
a form of finite A^ effects, if it is true, is in remarkable contrast to the "bosonic" case, in 
which the fermions in the IIB matrix model ( p.l| ) is omitted by hand. In that case, the large 
A^ behavior of correlation functions was determined analytically [|1^ to all orders in the 1/D 
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Figure 3: The 10 eigenvalues of the moment of inertia tensor normahzed by A^^/^ are plotted 
against IjN for N = 192, 256, 384, 512. 
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expansion and finite effects were found to be given by a l/N"^ expansion. 

For D = 6, the linear extrapolation gives almost the same results up to the third smallest 
eigenvalue, but not for the larger ones. Thus, as far as reproducing an isotropic distribution 
of Xi in the large limit is concerned, D = Q seems to have larger finite N effects. This 
rather counter-intuitive result might be understood by considering the branched-polymer 
representation of the low-energy effective theory The attractive potential between two 
Xj's connected by a bond is given by r~^^^~'^\ where r = a/ {xi — XjY- Therefore, it is 
much stronger for D = 10 than for D = 6. Since the branched-polymer system in both 
cases is not a simple maximal tree but something more complicated, typically involving a 
double-tree structure, it may well be that the stronger the attractive potential is, the more 
the distribution of Xi tends to become isotropic. 



5 Summary and Discussion 

In this paper, we studied the IIB matrix model (and its six-dimensional version) at large by 



using the low-energy effective theory developed in Ref. PSI. Unlike in the four- dimensional 



version studied in Ref. |T8[, the fermion integral yields a Pfaffian (or a determinant in the 
6D case) which is complex in general. In the present paper we omitted the phase of the 
Pfaffian by hand in order to avoid the sign problem in the Monte Carlo simulations, and we 
studied the effect of the modulus only. We have seen that the distribution of Xi becomes 
more and more isotropic as we increase N. Based on this observation we conclude that if 
the SSB ever occurs in the original IIB matrix model, the phase of the Pfaffian induced by 
fermions must play a crucial role. 

It is interesting to compare with the situation in perturbative superstring theory. Also 
in that case the fermionic degrees of freedom can be integrated out explicitly , leaving us 
with an effective bosonic string theory, where the action consists of three parts: the ordinary 
bosonic string action, an extrinsic curvature term, and (in Euclidean space-time) a purely 
imaginary Wess-Zumino-like term. This means that if we can neglect the imaginary part, the 
superstring theory is equivalent to a bosonic string theory with extrinsic curvature. While 
this equivalence has not been directly disproved, it does not have much support, either, from 



numerical simulations and analytical calculations [^, |29 



In Ref. [30 1 a deformation of the IIB matrix model by introducing an integer parameter 
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V which couples to the phase of the Pfaffian has been considered. The original IIB matrix 
model corresponds to z/ = 1. The deformed model continues to be well-defined, and preserves 
Lorentz invariance, the SU(A^) symmetry, and the cluster property. In this language the 
present work corresponds to the study of the case z/ = using the low-energy effective theory. 
The opposite extreme, the limit v = oo, has been studied analytically in Ref. pO|, and the 
spontaneous breakdown of Lorentz invariance has been discovered. Of particular interest 
is the fact that different conclusions have been obtained for u = oo and for u = 0, which 
already implies that a phase transition should occur in between. If the original IIB matrix 
model (z/ = 1) belongs to the same phase as z^ = oo, we can investigate the dynamics of the 
IIB matrix model by studying the z/ = oo model using Monte Carlo simulation as suggested 



in Ref. I^Oj- It would be extremely interesting if we can obtain a flat four-dimensional space- 
time in that way. We hope that the Monte Carlo technique developed in the present work is 
useful for such future studies and eventually enables us to explore the dynamics of the IIB 
matrix model. 
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Appendix A: The algorithm for the simulation 



In this appendix, we explain the algorithm]^ we use for the Monte Carlo simulation of our 
model defined by { ^2^ . 

We first simplify the model by introducing rescaled variables fo^jj as 



The model we have to simulate then becomes 



where 



Z = y"dxd6e-^'3-^^"' |PfM| , (A.2) 



SG = Y.\^,^^\ (A.3) 



and ( |3.12| ) can be rewritten in terms of hij as 



ML,, = ^(^^ l"";^.), ^^(|^)2|2 (r/^r.).,(&,,,6.., - for z^j, (A.4) 

where we have also introduced the dimensionless regularization parameter e to reduce nu- 
merical instabilities. The regularization parameter e should be taken to be small enough 
not to affect the system. For both, 6D and lOD, cases, we took e = 0.1. We note also 
that with this regularization, the model is well-defined without introducing the cutoff term 
( p.l9| ). However, by simulating such a model, we find that all the Xj's become densely packed 



within the extent e^ as we increase A^. Therefore, we still need the cutoff term ( p.l9| ) in order 
to make our model valid as a low-energy effective theory of the original IIB matrix model. 
For convenience we combine our dynamical variables Xi^ and {i 7^ j, = 6*^^) into a 
Hermitian matrix as {B^)ij = Xi^5ij + h^ij. 

The algorithm we use to simulate the model ( [A. 21 ) is a variant of the Hybrid Monte Carlo 



algorithm ||2^. The first step of the Hybrid Monte Carlo algorithm is to apply molecular 
dynamics [Q. We introduce a conjugate momentum for B^ij as X^ij, which satisfies X^ji = 
^Ref. gives an overview of effective algorithms for dynamical fermions. 
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(X/j^ij)*; are Hermitian matrices. The partition function can be rewritten asQ 

Z = j dXd5e-^ , (A.5) 
where H is the "Hamiltonian" defined by 

H=\Y. I^/^^^-I' + "^gI^] + 5cut[5] - i In I det M| . (A.6) 
m 

The update of X^jj can be done by just generating X^^- with the probabihty distribution 
exp(— i Yli l-'^wH- -'^^ order to update -B^jj, we solve the Hamilton equation given by 



dR„.fr) dH 



dr dX^ij 



X,,^ (A.7) 



dX,.,(r) ^ ^ iTr r^D-i V ^ - ^ fA 8) 

dr 95^,, 4 V55^,, ) dB^,, dB^,, ' ^ ■ ^ 

where we have defined D = M'^ M. Note that when taking the derivatives in ( |A.8| ), {B^)ij and 
{B^)ji should be treated as independent complex variables. Along the "classical trajectory" 
given by the Hamilton equation, 

(i) H is invariant, 

(ii) the motion is reversible, 

(iii) the phase- volume is preserved; i.e. 

3(B(r),X(r)) ^ (A 9) 

S(i3(0).A-(0)) ' 

where {B{t),X{t)) is the point on the trajectory after evolution for fixed r from {B{0), X{0)). 
Therefore, generating a new sets of {B,X) by solving the Hamilton equation for a fixed 
"time" interval r satisfies the detailed balance. This procedure, together with the pro- 
ceeding generation of X^ij with the Gaussian distribution, is called "one trajectory", which 
corresponds to "one sweep" in ordinary Monte Carlo simulations. 

When solving the Hamilton equation numerically, we have to discretize the equation. A 
discretization which preserves the properties (ii) and (iii) is known. The property (i) cannot 
be preserved and yields a small violation of the Hamiltonian conservation. In order to still 
^Unlike in the standard Hybrid Monte Carlo algorithm p6|, we do not introduce the so-called pseudo- 



fermions. Note, in this regard, that our system (A. 2) is different from ordinary field theories with dynamical 
fermions in the following respects. 1) The number of fermion flavors should be strictly one in order to respect 
supersymmetry. 2) The size of the matrix M is much smaller than the system size. 3) The matrix M is 
not sparse. 
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satisfy the detailed balance exactly, we can perform a Metropolis accept/reject as the end 
of each trajectory. 

We introduce a short-hand notation for discretized X^{t) and B^ij) as 

XM = X,(rAr) ; 5^ = i?,(sAr) . (A.IO) 

Given the configuration {Bf^\x^'^), the configuration {B^^\x\^'^) after the evolution for a 
fixed time r = v At is obtained by solving the discretized Hamilton equation 

bI:;^' = B^'^+^rx^ 

= 4S-ArJ^(5r^^), (A.12) 



where n = 1, 2, ■ ■ ■ , z/— 1 and m = 0, 1, ■ ■ ■ , z/— 1. Note that the first and the final steps in the 
evolution ( A. 11 ) of -B^ are treated with special care. This particular discretization, which 
is called as "leap-frog discretization", preserves the properties (ii) and (iii). At the end of 
the trajectory, we make a Metropolis accept/reject with the probability P = min(l, e"'^''^), 
where AH is the difference of the Hamiltonian H defined in ( |A.6| ) for the configurations 
{Bj^\xi^^) and {B^^\xIj:^). Optimization of Ar and V is discussed in Appendix B. 
In the evolution ( |A.12| ) of X^, one needs to calculate 

If we write the first term as T^jj, the second term can be written as T^ji- In what follows, 
we calculate T^ij explicitly. We first note that 



k,l=l al3 f^""^ k,l=l al3 ^""^ 

where Cia^jp is defined as 

Cic,,jp = (M-^)i,j^ , (A.15) 

and Cl^ j^ is defined as 

Ct /^i \ ^ 

ia,jl3 — ^iaJP j ^ia,Np — ~ ^ia,k^ 

k=l 

N-1 N-l 

^NajP = ^ ^ Cka,jl3 ] C'j^a,N/3 = ^ Cka,ll3 ■ (A. 16) 

k=l kJ=l 



r)M' 

k=l l^k a/3 ™ 



The indices i,j in ( |A.15| ) and ( |A.16D run from 1 to — 1. We further rewrite T^ij as 

N N 

dB 

where we have introduced C'l^ ^f^ through 

Then we calculate ^^"''^ explicitly, which yields 



- _ ^ y _^^^g£-)2j2'f^^'^^l'^/^)°/5 i^p.^l^'T)al3}KjiC"^ j^ (A.19) 

for 2 7^ J and 

(r<TrprT)a/3(&o-fci&rifc — b„ikbrki)C"^ka ; (A. 20) 

where we defined Clajp through 

= ^ia,jP ~ ^ja,jl3 " Cia,if3 + ^ja,i0 ■ (A.21) 

Let us comment on the required computational effort of our algorithm. The dominant 
part comes from calculating the inverse in ( |A.15| ), which requires O(n^) arithmetic opera- 
tions, where n is the size of the matrix to be inverted. In the present case n is of O(A^). In 
order to keep the acceptance rate at the Metropolis accept/reject procedure reasonably high, 
one has to decrease the step size At as one goes to larger A^. We have seen from simulations 
that the step size should be taken to be At ~ which is consistent with a general formula 
At ~ V~^^^ in Ref. ||3^, where the system size V should be replaced by 0(A^^) in our case. 
Accordingly, the number of steps for one trajectory increases as z/ ~ Therefore, the 

required computational effort is estimated to be 0{N'^^'^). 

In fact, one can reduce the computational effort by omitting the Metropolis accept/reject 
procedure, since in that case one can keep the step size At fixed to a small constant for all 
A^. The required computational effort becomes 0{N^). The price one has to pay is that the 
algorithm then suffers from a systematic error due to the small violation of the Hamiltonian 
conservation. The systematic error can be estimated to be 0((Ar)^), see Appendix B. 
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Appendix B: Optimization of the algorithm 



In this appendix, we first discuss the optimization of the parameters in the algorithm with 
the Metropohs accept/reject procedure. We then move on to the case when the Metropohs 
accept /reject procedure is omitted. The algorithm and the parameters used for each run is 
also described. 

To start with, one can actually generalize the algorithm described in Appendix A by 
taking the step size At^ for the diagonal elements to be different^ from the step size Ar for 
the off-diagonal elements of B^ij and X^ij. We have fixed the ratio Ar^;/ Ar to be the ratio of 
the mean magnitudes of Xi^ and which we found from simulations to be approximately 
1 for D = Q and 1/2 for D = 10, respectively. 

We still have two parameters in the algorithm, i.e., Ar and v (the number of molecular 
dynamics steps for each trajectory), which can be optimized in a standard way [35, 36|. The 



key point of the optimization is that the autocorrelation time (in units of accepted trajectory) 
depends only on v Ar = r, but not on Ar and v separately. This allows us to perform the 
optimization in two steps. First, one fixes v At = r and optimize Ar so that the effective 
speed of motion in the configuration space, given by acceptance rate times Ar, is maximized. 
Second, using the optimized Ar for each r, one minimizes a typical autocorrelation time in 
units of molecular dynamics step with respect to r. 

In Fig. ^, we show the acceptance rate times Ar as a function of Ar ioT D = 6, N = 16 
with fixed r = 1.5. We find that the optimal Ar is 0.0375. For D = 6,N = 32 with 
fixed r = 1.5, we find that the optimal Ar is 0.02142, which is smaller than for = 16 as 
expected. For both cases, the acceptance rate at the optimal Ar is found to be 50 ~ 60%. 

Using the optimal Ar obtained in the above way for each r, we minimize a typical 

autocorrelation time (in units of molecular dynamics step) with respect to r. Here, we 

measure the autocorrelation time of the extent R of the Xj-distribution defined in ( |4.2D and 

plot it as a function of r. Fig. ^ shows the result for D = 6 with A^ = 16. We see that it 

has a minimum around r ~ 1.5. Similar experiments for A^ = 32 showed that the optimal r 

is almost independent of A^. 

®This is not the case for the fuU model where At should be taken universally for each element of 
A^ij in order to respect the SU(iV) invariance. In the present case, since the SU(A^) "gauge" invariance is 



fixed by the "gauge-fixing" (3.4), we only have to respect the invariance under permutation of the SU(A^) 
indices. 
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Figure 5: The acceptance rate times At is plotted against At ior D — 6, N — 16 with fixed 
T = 1.5. 
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When we omit the Metropohs accept/reject procedure, the algorithm suffers from a 
systematic error, as is explained in Appendix A. Here we fix v Ar to the optimal r obtained 
for the case with the Metropolis accept/reject procedure and study the Ar dependence of 
the systematic error. As a quantity which shows a large systematic error, we take (Ai), the 
expectation value of the largest eigenvalue of the moment of inertia tensor defined by ( ^.1| ). 
In Fig. 1^ we plot the result against (Ar)^ for D = 6, = 16 with r = 1.5. The systematic 
error is seen to vanish as 0((Ar)^), which can be also understood theoretically by using the 
analysis described in Ref. |l36| . 

4 




Figure 7: The expectation value of the largest eigenvalue of the moment of inertia tensor is 
measured by simulations without Metropolis accept/reject procedure. The result is plotted 
against (Ar)^ for D = 6, N = 16 with fixed r = 1.5. The filled circle at Ar = represents 
the result obtained by a simulation with Metropolis accept/reject procedure. The straight 
line is a fit to the predicted (Ar)^ behavior of the systematic error. 



Finally, let us comment on the algorithms and the parameters we used in our simulations 
at large A^. The runs for D = 6 with A^ = 192, 256 were made by the algorithm including 
the Metropolis accept/reject procedure. The parameters are Ar = 0.0048, = 60 for A^ = 
192, and Ar = 0.004, z/ = 50 for A^ = 256. The numbers of configurations used for the 
measurements are 560 and 1732 for A^ = 192 and A^ = 256, respectively. (The optimization 
described in this Appendix was not completed when we started these runs. Accordingly, we 
needed significantly larger numbers of configurations compared to the cases below.) For the 
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other cases, we omitted the Metropohs accept/reject procedure in order to obtain a sufficient 
statistics. The parameters are Ar = 0.0075, v = 200 for D = 10 with = 192, 256, 384, 512, 
and Ar = 0.005, u = 300 for D = 6 with N = 512, 768. The numbers of configurations used 
for the measurements are 88, 40, 84, 44 for D = 10 with = 192, 256, 384, 512, and 64, 50 
for D = 6 with N = 512, 768, respectively. 
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